Wave functions for quantum Monte Carlo calculations in solids: Orbitals from density 
functional theory with hybrid exchange-correlation functionals 
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We investigate how the fixed-node diffusion Monte Carlo energy of solids depends on single- 
particle orbitals used in Slater-Jastrow wave functions. We demonstrate that the dependence can be 
significant, in particular in the case of 3d transition-metal compounds, which we adopt as examples. 
We illustrate how exchange-correlation functionals with variable exact-exchange component can be 
exploited to reduce the fixed-node errors. On the basis of these results we argue that the fixed-node 
quantum Monte Carlo provides a variational approach for optimization of effective hamiltonians 
with parameters. 
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I. INTRODUCTION 

Steady increase in supercomputer performance over 
the last three decades stimulates development of very ac- 
curate electronic structure methodologies that provide 
detailed understanding of ever larger and more com- 
plex systems. Diffusion Monte Carlo (DMC) method 1,2 
is one of the prime examples of these advanced many- 
body approaches. It uses a stochastic process to re- 
fine a given many-body wave function \^t) towards the 
actual ground-state solution |^o)- The antisymmetry 
of the wave function with respect to fermionic degrees 
of freedom is usually maintained by imposing the so- 
called fixed-node condition. In practice it means that the 
nodal surface (the subset of configuration space where the 
many-body wave function vanishes) is restricted to be the 
same as in the initial guess \^t) throughout the entire 
stochastic simulation. The basic premise is that high ac- 
curacy is achieved even when relatively simple functional 
forms are employed for \^t)- 

The impact of the fixed-node approximation has been 
intensively studied in few-electron atoms and molecules 
as well as in homogeneous systems, where the perfor- 
mance of the DMC method with increasing accuracy of 
the trial wave function \*&t) is relatively well mapped 
outj^~— Applications to solids are not nearly as numer- 
ous, and hence the influence of the fixed-node condition 
in crystalline settings is much less examined. Evalu- 
ation of bulk properties necessarily involves extrapola- 
tion to the thermodynamic limit, which decreases the 
amount of computational resources available for explo- 
ration of rather subtle fixed-node errors. Consequently, 
only the simplest trial wave functions, having the Slater- 
Jastrow form, are typically employed. These wave func- 
tions have nodal surfaces defined by single-particle or- 
bitals, and therefore they correspond to mean-field nodes. 
The subject of our study are the fixed-node errors asso- 
ciated with such trial functions when they are applied to 
simple transition-metal oxides in crystalline phases. 



II. FIXED-NODE DIFFUSION MONTE CARLO 

The fixed-node diffusion Monte Carlo method calcu- 
lates expectation values of quantum-mechanical opera- 
tors in the ground state (as well as in certain excited 
states) of a many-body hamiltonian H . The DMC wave 
function is found by a projection 
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that gradually increases the weight of the lowest energy 
eigenstate of H relative to all other states contained 
in some initially guessed wave function \^t)- Build- 
ing on similarity between the Schrodinger and the dif- 
fusion equations, the projection is realized with the aid 
of a classical stochastic process. The outcome of this 
simulation is a set of 3iV-dimensional samples {Ri} dis- 
tributed according to a probability distribution V(R) = 
*d(R)*t(.R)/{*d\*t). Here R = (n,ra, . . . , r N ) de- 
notes coordinates of all N electrons comprising the inves- 
tigated quantum system. The hamiltonian H is assumed 
to be spin independent, which prevents any spin- flip pro- 
cesses to occur in Eq. (fTa)) . and hence the spins do not 
enter the simulations as dynamical variables. The prob- 
abilistic interpretation of V(R) is possible only if it is a 
positive quantity. In the case of fermions, the projection 
as written in Eq. (|la[) does not fulfill this requirement, 
which must therefore be prescribed in the form of an ad- 
ditional condition 



y D (R)y T (R) > 0. 
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This step introduces the so-called fixed-node approxima- 
tion, since the projection cannot reach the true ground 
state |^o) if the fermionic nodes of the trial wave func- 
tion I^t) differ from (a priori unknown) nodes of the 
desired ground state. Our primary aim is to explore the 
impact of this fixed-node approximation for a particular 
functional form of the trial wave function. 

The set of J\f samples {Ri} acquired as a result of a 
DMC simulation can be directly used to calculate the 



2 



so-called mixed estimates of the quantum-mechanical ex- 
pectation values, 
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If the operator A commutes with the hamiltonian, 
the mixed estimate equals the desired quantum- 
mechanical expectation value, {^d\A\^>t)/{^ d\^t) = 
d\A\^ d) / i 1 ^ d\^ d) ■ In general, however, there is an 
error proportional to the difference between \^d) and 
\^>t) that can be reduced to the next order using the 
following extrapolation 2 - 
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The expectation value (^tIA^t) / (^tI'&t) is evaluated 
by straightforward Monte Carlo integration and it is re- 
ferred to as the estimate of the variational Monte Carlo 
(VMC) method. 



III. TRIAL WAVE FUNCTIONS 

We employ trial wave functions having the Slater- 
Jastrow functional form, which is an antisymmetrized 
product of single-particle orbitals (Slater determinant) 
multiplied by a correlation factor that is symmetric with 
respect to pair-electron exchanges. We can write 



V T (R) = det{^} det {^} e J ^ x > , 



(4) 



where and ijA are spatial parts of single-particle or- 
bitals that correspond to spin-up respectively spin-down 
electronic states. The vector X = (xi,X2, . . . ,Xm) en- 
compasses positions of all M ions in the lattice. The 
expression in Eq. Q represents only a spatial compo- 
nent of the trial wave function corresponding to one par- 
ticular spin configuration, where electrons with labels 
1, . . . , iV-f are in the spin- up state and electrons with la- 
bels iV-f + 1, . . . , N are in the spin-down state. We can use 
this simplified form with fixed spin states instead of the 
full wave function as long as neither the hamiltonian nor 
any other operator, expectation value of which we intend 
to calculate, depend on electron spins. The applications 
we consider involve only cases with zero total spin, i.e., 
N is an even number and Nf — N± = N/2. 
The Jastrow correlation factor we use, 

J{R, X)=J2 f(*> - rj ) + E ^ - x ') ' ( 5 ) 



contains one- and two-body terms, g and /, that are 
parametrized in the same way as in Ref. H|. This cor- 
relation factor improves efficiency of the Monte Carlo 
sampling and accuracy of general expectation values cal- 
culated according to Eq. ([3|). Quality of the DMC total 
energy Ed depends solely on the accuracy of the nodal 
surface that is, given the functional form of Eq. (TJ}, 
completely determined by the single-particle orbitals 
V'jj}- Ideally, these orbitals would be parametrized 
by an expansion in a saturated basis with the expan- 
sion coefficients varied to minimize the DMC total en- 
ergy. Unfortunately, the stochastic noise and the compu- 
tational demands of the DMC method make this route 
extremely inefficient in practice except, perhaps, in the 
case of the simplest few electron systems. 

A more feasible alternative is to skip the DMC pro- 
jection, Eq. (jlal) . and optimize the orbitals {VaiV^i 
with respect to a simpler quantity than Ed- For in- 
stance, the variational Monte Carlo methodology can 
be employed to minimize the variational energy Ey = 
{^ t \H\^ t )/(^t\^t}- The VMC optimization of one- 
particle orbitals was successfully employed for atoms and 
small molecules of the first-row atoms ^ 10 ' 11 but the 
method is still too demanding for applications to solids. 
It is also not completely robust, since an improvement of 
the variational energy does not automatically guarantee 
an improvement of the fixed-node energy due to the lim- 
ited parametric freedom of the trial wave function \^t)- 
This issue is even more pronounced when trial wave func- 
tions are not optimized with respect to the VMC energy, 
but with respect to another quantity, such as the energy 
variance. 

To avoid the large number of variational parameters 
needed to describe the single-particle orbitals, another 
family of methods has been proposed. The orbitals in the 
Slater- Jastrow wave function are found as solutions to 
self-consistent-ficld (SCF) equations that represent a gen- 
eralization of the Hartree-Fock equations! 12 ! 13 The cor- 
relations described by the Jastrow factor enter the SCF 
equations via variational Monte Carlo techniques. These 
methods were tested in atoms as well as in solidS ) 12 i 14 al- 
beit only rather limited variational freedom was allowed 
in the employed Jastrow factors. Unfortunately, wave 
functions derived in this way did not lead to lower DMC 
energies compared to wave functions with orbitals from 
the Hartree-Fock theory or from the local density ap- 
proximation (LDA) . 12 ' 15 It is unclear, whether the lack 
of observed improvements in the fermionic nodal surfaces 
stems from insufficient flexibility of the employed corre- 
lation factors or from too small set of tested cases, all 
with only s and p valence electrons. 

In this paper we also use SCF equations as a means 
to construct the one-particle orbitals, but the parametric 
dependence of these equations is introduced without any 
relation to the actual Jastrow factor. The self-consistent- 
field equations in question are Kohn-Sham equations cor- 
responding to the hybrid exchange-correlation functional 
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PBEO^ 

£™» = wE BF + (1 - w)E PBE + E P C BE . (6) 

Here E PBE and E PBE are exchange and correlation 
parts of the PBE formal of the generalized gradient 
approximation (GGA), and E BF is the exchange from 
the Hartree-Fock theory. The weight w is in the range 
< w < 1 and serves as a variational parameter with re- 
spect to which the fixed-node DMC energy is minimized. 

In the transition-metal oxides we study it is understood 
that exchange in GGA is underestimated^ whereas the 
"exact" exchange from the Hartrcc-Fock theory overes- 
timates the real exchange mechanism as any screening 
effects are neglected i 19 ' 20 The hybrid density-functional 
theory (DFT) provides an interpolating scheme between 
these two extremes. There are numerous examples in 
the literature illustrating that reasonable agreement of 
its predictions and experimental observations can be 
achieved £±r— 

Arguably, the single parameter introduced to the one- 
particle orbitals represents only relatively constrained 
variational freedom compared to the approaches men- 
tioned abovei 6 ' 10 ' 11 On the other hand, the simplicity of 
the parameter space allows for direct optimization of the 
fixed-node DMC total energy, which improves robustness 
compared to minimization within the VMC method. Ini- 
tial applications of this strategy to molecules 9 -^ and in a 
more elementary form also to solids^ 6 - were reported pre- 
viously. In the following sections we perform the DMC 
optimization for two compounds, MnO and FeO, and 
systematically analyze the variational freedom available 
within this method in crystalline environments. 

IV. THERMODYNAMIC LIMIT 

Crystals in our simulations are represented by period- 
ically repeated simulation cells of a finite size, in which 
the Coulomb interaction energy is evaluated with the aid 
of the Ewald formulai^ To calculate bulk quantities one 
has to perform an extrapolation to the thermodynamic 
limit (infinite crystal volume). Our main objective, how- 
ever, is to compare energies obtained for different trial 
wave functions in a given simulation cell with fixed vol- 
ume, and hence the size extrapolation seems redundant. 
Indeed, the leading-order term of the finite-size errors re- 
lated to the long-range character of the Coulomb interac- 
tion is a function of the average charge density alonei 28 ' 29 
Variation of single-particle orbitals does not change this 
quantity, and these finite-size errors are therefore mostly 
irrelevant for our wave function optimization, since they 
cancel out. 

The Coulomb interaction is not the only cause of finite- 
size biases in calculations performed in finite simulation 
cells. There is another source of non- vanishing surface 
terms that appear even if no two-body interaction is 
present in the hamiltonian. 30 These can be brought to 
light as follows: the statement "periodically repeated 



TABLE I: Exact-exchange weight w leading to minimal DMC 
total energies of several phases of MnO and FeO solids. 



compound phase 


v (A 3 /xo) 


w 


MnO 


Bl AFM-II 


21.7 


0.30 ±0.01 


MnO 


B8 AFM 


21.7 


0.33 ± 0.01 


FeO 


Bl AFM-II 


20.4 


0.25 ± 0.02 


FeO 


Bl AFM-II 


17.3 


0.217 ±0.006 


FeO 


iB8 AFM 


17.0 


0.156 ±0.003 



simulation cells" used above refers to observable quan- 
tities and does not fully specify the boundary condi- 
tion for the phase of the wave function. In the lan- 
guage of single-particle description of solids, different 
boundary conditions compatible with periodicity of ob- 
servables correspond to simulations being done at dif- 
ferent k points. Averaging over all k points from the 
first Brillouin zone removes the dependence of the com- 
puted quantities on the boundary conditions and in the 
case of non-interacting particles it is equivalent to per- 
forming the thermodynamic limit. For hamiltonians with 
particle-particle interactions this correspondence is not 
exact, but a very substantial reduction of finite-size er- 
rors is observed nevertheless^ - 

Comparison of trial wave functions at a single k point is 
certainly a valid approach. On the other hand, variation 
of the exchange-correlation functional in the Kohn-Sham 
equations modifies the resulting band structure in a gen- 
erally non-trivial fc-dependent manner, see for instance 
Ref. l2l| for illustration. As a result, the minimization of 
the total energy calculated at a single k point does not 
necessarily lead to the same optimal weight w as does the 
minimization of the fc-averaged total energy. Since the 
latter is arguably a better approximation of the thermo- 
dynamic limit, which is what we are ultimately interested 
in, we work with the fc-averaged quantities. 



V. MnO AT AMBIENT CONDITIONS 

We start our investigations with manganese monoxide 
at experimental equilibrium volume, V = 21.7 A 3 /MnO. 
Calculations are performed for two lattice structures, 
both with antiferromagnetic (AFM) ordering of mag- 
netic Mn atoms: Bl (symmetry group R3m, the so-called 
AFM-II state) and B8 (symmetry group P3ml). The for- 
mer phase is the low-temperature ground state structure 
at atmospheric pressure and the latter is a phase stable 
at high pressures^ Since all states investigated in this 
paper are antiferromagnetic, we will often leave out this 
attribute. 

Data were collected in simulation cells containing 16 
atoms (8 Mn and 8 O). Atomic cores were replaced by 
norm-conserving pseudopotentials within the so-called 
localization approximation. 32 Helium core was excluded 
from oxygen atoms^ 3 - and neon core from manganese 
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FIG. 1: (color online) The DMC total energy of MnO at 
experimental equilibrium volume (V = 21.7 A 3 /MnO) for 
two distinct structural phases: Bl AFM-II (red squares, solid 
line) and B8 AFM (blue points, dashed line). Minima of the 
fits are chosen as reference energies for the respective phases. 
Also shown is the DMC total energy of the Bl AFM-II phase 
obtained with the Hartree-Fock orbitals (empty red square). 



FIG. 2: (color online) The DMC total energy of MnO (Bl 
AFM-II phase at V = 21.7 A 3 /MnO) as a function of exact- 
exchange weight w. Compared are data obtained in simula- 
tion cells of two different sizes: 16 atoms (red squares, solid 
line) and 4 atoms (blue circles, dashed line). The lines are 
quadratic fits. Minima of these fits are chosen as reference 
energies for the respective data sets. 



atoms; 34 which left 168 valence and semi-core electrons 
explicitly included in the simulations. 

The Monte Carlo simulations were performed with 
QWalk££. and the single-particle orbitals were prepared 
in Crystal2003i2£ In both codes, the orbitals were ex- 
panded in a gaussian basis, completeness of which was 
verified against a converged basis of linearized augmented 
plane waves as implemented in the Wien2k code^i The 
basis set benchmarks were performed within DFT using 
the PBE-GGA exchange-correlation functional. 

The acquired dependence of the DMC total energies on 
the exact-exchange weight w is plotted in Fig. [TJ Each 
energy value shown in the picture represents an average 
over 8 k points. Such average made within DFT dif- 
fers from the converged integral over the first Brillouin 
zone by less than the statistical error bars plotted in the 
figure. The quadratic functions fitted through the data 
using the least-squares method lead to the optimal values 
for the weight w listed in Tab.Q] The behavior is virtually 
the same in both phases and the optimal weight is very 
close to the standard PBE0 choice^ w = 0.25. Drop 
of the total energy from the pure PBE-GGA orbitals 
(<=> PBE0„ = o) to the minimum is w 0.15 eV/MnO, and 
from the Hartree-Fock orbitals (<ft PBE0 w=1 ) it is w 0.2 
eV/MnO. The dependence of the fixed-node DMC total 
energy on the admixture of the exact exchange is weaker 
in the MnO solid than in the MnO molecule, where the 
gain from the Hartree-Fock to hybrid DFT orbitals was 
observed as large as ~ 0.6 eV<2 We attribute this differ- 



ence to a greater freedom of charge density to adapt in 
the molecule than in the constrained solid state environ- 
ment. 

The data collected so far can be utilized to evaluate 
quantities of direct physical interest. The difference be- 
tween the optimized total energies provides a measure of 
relative stability of the investigated phases at the given 
volume. We find that the Bl structure is lower in en- 
ergy than the B8 structure by 0.27 ±0.01 eV/MnO. The 
cohesive energy E co h of MnO crystal can also be esti- 
mated. To this end, energies of isolated Mn and O atoms 
need to be calculated, and the crystalline total energy 
has to be extrapolated to the infinite volume. We pro- 
vide details of these steps in Appendix [3] In the end we 
obtain E co h — 9.29 ± 0.04 eV that is in good agreement 
with the value 9.5 eV derived from experimental forma- 
tion enthalpies.— Our result differs from the earlier DMC 
estimate ; 39 ' 40 which is not surprising as the present ap- 
proach involves more advanced optimization of the trial 
wave function as well as much more thorough analysis of 
the finite-size errors. 



VI. ARE LARGE SIMULATION CELLS 
NECESSARY? 

It was argued in Sec. HVI that k point averaging should 
remove nearly all finite-size biases relevant to the opti- 
mization of single-particle orbitals performed at a fixed 
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volume. It would certainly be beneficial if one could find 
the optimal exact-exchange weight in a small simulation 
cell and only then proceed with production runs in large 
cells. 

We have repeated the optimization procedure in the 
primitive cell of the Bl AFM-II structure, which contains 
only 4 atoms (2 Mn and 2 0). All other parameters were 
kept unchanged, except the number of k points employed 
in the averaging had to be substantially increased — from 
8 to 125 — to achieve comparable convergence.— The ne- 
cessity to enlarge the set of considered k points intro- 
duces an extra technical complication, since trial wave 
functions corresponding to the majority of these k points 
are complex-valued. Consequently, we replace the fixed- 
node DMC method with its natural generalization, the 
so-called fixed-phase DMC^i The fixed-phase condition 
reduces to the fixed-node condition for real-valued trial 
wave functions. 

The w-dependence of the DMC total energy calculated 
in the small simulation cell (Fig. [5]) is similar but not com- 
pletely identical to the behavior observed in the larger 
cell. The minimum is shifted to a slightly lower weight w 
and the E(w) curve raises a little slower with w decreas- 
ing towards the pure PBE-GGA. It is plausible to assume 
that these differences are a fingerprint of the residual 
Coulomb finite-size effects. 



VII. COMPRESSED FeO 

Finally, we turn our attention to iron oxide subject to 
high pressures at low temperatures. We revisit our earlier 
study of the phase transition from the Bl AFM-II phase, 
stable at atmospheric pressure, to the iB8 AFM structure 
(space group P6m2), which the fixed-node DMC has pre- 
dicted as stable above approximately 65 GPa. 42 Calcula- 
tions leading to this estimate used single-particle orbitals 
provided by the PBEO„, = o.2 functional. This choice was 
based on prior investigation !] 9 ' 25 and on a preliminary 
version of the results we present here. In the following, 
we analyze in detail how appropriate the PBEOu, = o.2 or- 
bitals are in this case and how sensitive the transition 
pressure is to variations of single-particle orbitals in the 
Slater- Jastrow trial wave function. 

We return back to simulation cells containing 16 atoms, 
partly because the orbital optimization appears more ro- 
bust in larger cells (Sec. IVI|) and partly because we al- 
ready acquired some data in the 16 atom cells in the 
course of our earlier investigations. The optimized val- 
ues for the exact-exchange weight w in FeO are listed 
in Tab. U for experimental equilibrium volume V = 20.4 
A 3 /FeO as well as for two compressed states. The opti- 
mal proportion of the exact exchange decreases with com- 
pression, which is an expected phenomenon — the role of 
screening increases as the bands widen and a larger frac- 
tion of d electrons participates in chemical bonding. 

In contrast to MnO, where Bl and B8 phases displayed 
very similar behavior, the two investigated structures of 



0.30 




-0.05 1 1 1 1 1 1 

0.0 0.1 0.2 0.3 0.4 0.5 
exact exchange weight in PBE0 w 



FIG. 3: (color online) The DMC total energy of compressed 
FeO in two phases: Bl AFM-II at V = 17.3 A 3 /FeO (red 
squares), and iB8 AFM at V = 17.0 A 3 /FeO (blue points). 
Lines are least square fits with a quadratic function in the 
iB8 phase, and with E(w) = (a + bw + cw 2 )/(d + w) in the 
Bl structure. Note the different choice of reference energies 
compared to Figs. [T] and O Here the zero energy for each 
phase is the value of the corresponding fit at w = 0.2. 



FeO differ noticeably at comparable volumes. Detailed 
data for compressed FeO are shown in Fig. [3] to high- 
light the differences. Like in Sec. [V] each total energy is 
obtained as an average over 8 k points. To facilitate com- 
parison with our previous study ^2 energies correspond- 
ing to the PBEO^o.2 orbitals are used as a reference. 
We can see that w = 0.2 indeed represents a reasonable 
compromise value, since corresponding DMC energies lie 
within error bars from the true minima. 

In the case of the iB8 phase, the E(w) data are well 
described by a quadratic function. In the Bl phase, 
on the other hand, the functional dependence is asym- 
metric around the minimum and a quadratic function 
does not provide a satisfactory fit. To locate the mini- 
mum for Tab. [J we used an alternative fitting function 
E(w) = (a + bw + cw 2 )/(d + w), which characterizes the 
calculated energies much better. Figure [3] shows that the 
energy raises rather fast when the exchange-correlation 
functional approaches the pure PBE-GGA. It does not 
come as a surprise, since the Kohn-Sham spectrum is 
metallic in the limit w — > 0, which is at odds with exper- 
imental facts. All other phases investigated here (FeO 
iB8 and both MnO phases) are insulating for all values 
of the weight w, and hence even orbitals close to the 
PBE-GGA provide reasonable trial wave functions. 

The different behavior of the DMC total energies in the 
Bl and iB8 structures causes the pressure of the transi- 
tion between these two phases to depend on the used 
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orbitals. To roughly estimate the variation of the transi- 
tion pressure, we assume that the energy- volume equa- 
tions of state of the respective phases only rigidly shift 
along the energy axis when the exchange weight is varied. 
The equations of state corresponding to the PBE0u, = o.2 
orbitals were calculated in Ref. and the correspond- 
ing shifts can be extracted from Fig. [3J When, for in- 
stance, the PBEO^o.3 orbitals are used, the iB8 phase 
is raised in energy by approximately 0.1 eV compared to 
the Bl structure, which leads to the transition pressure 
increased to as 85 GPa. When the PBEO^^o.05 orbitals 
are utilized, the iB8 phase is lowered by approximately 
0.2 eV, which corresponds to the transition pressure of 
only as 30 GPa. Evidently, the Bl to iB8 transition pres- 
sure is quite sensitive to the choice of the single-particle 
orbitals. Of course, the DMC method provides a definite 
prediction as long as the exchange weight is individually 
optimized in each phase. The pressure 65±5 GPa derived 
in Ref. [42J remains valid as the actual DMC estimate for 
the Bl to iB8 transition pressure, since the energies ob- 
tained with the PBEOm^o.2 orbitals lie within error bars 
from the minimal energies (Fig. 13]). 



VIII. OPTIMIZATION OF EFFECTIVE 
HAMILTONIANS 



The optimization of not only the variational wave func- 
tion but also of the effective hamiltonian can be consid- 
ered in a broader context. The upper-bound property 
of the fixed-node approximation offers a new opportu- 
nity for finding optimal values of any parameters of effec- 
tive hamiltonians in a consistent and well-defined man- 
ner. For example in DFT+U methods^ the Hubbard 
parameter U as well as the form of the double counting 
terms could be optimized in a similar way as the exact- 
exchange weight in the presented calculations. The opti- 
mized hamiltonian can subsequently be utilized for fur- 
ther calculations. We illustrate this on the density of 
states of the iron oxide shown as a function of the exact- 
exchange weight in Fig. 2) Note that for PBEO„, = o.2, he., 
for the functional that we used as optimal for the calcula- 
tions of the equation of state, the single-particle spectrum 
exhibits a reasonable value of the gap of about 2.9 eV, 
which is much closer to the experimental value (as 2.4 
eV, Ref. |43T ) than either pure GGA or pure Hartree-Fock 
limits. The key point is that although the hamiltonian 
was optimized in the ground state, the excitations are 
clearly improved as well. 

In a more general sense, one could consider effective 
hamiltonians with more parameters and/or with more 
elaborated content beyond the one-particle form, such 
as explicit treatment of particle pairs, for example. The 
idea of employing the fixed-node DMC method could be 
of importance for such constructions as the most accurate 
and explicitly variational method available at present. 
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FIG. 4: (color online) Density of states of FeO (Bl phase at 
V = 20.4 A 3 /FeO) from PBE0™ calculations with varying w. 
Where it is not obvious, the bottom of the conduction band 
is indicated with an arrow. 



IX. CONCLUSIONS 

We have found that single-particle orbitals in Slater- 
Jastrow wave functions represent a non-trivial variational 
parameters for fermionic nodes in 3d transition-metal 
compounds. When these orbitals are generated with 
the aid of an exchange-correlation functional with vari- 
able admixture of the exact exchange, the corresponding 
fixed-node DMC energies differ by several tenths of eV 
per transition-metal atom. These variations can trans- 
late to substantial quantitative changes in the phase di- 
agram as demonstrated in the case of iron oxide. Con- 
sequently, some form of orbital optimization should be 
performed in order to confidently predict relative sta- 
bility of different crystal structures and the location of 
corresponding phase transitions. 

The optimal amount of exact exchange providing min- 
imal DMC energies depends on the particular compound 
and on the specific structural phase, but it is generally 
close to the value of 25% typically used within the hy- 
brid density-functional theoryJ^ Our calculations pro- 
vide yet another confirmation that the hybrid DFT in- 
deed provides an improved picture of investigated com- 
pounds compared to more conventional functionals (LDA 
and GGA). 
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16 12 8 
No. of MnO formula units 

FIG. 5: (color online) Finite-size scaling for MnO (Bl AFM- 
II phase at V = 21.7 A 3 /MnO). All data points are averaged 
over 8 k points. Statistical error bars are smaller than symbol 
sizes. Red squares correspond to the Ewald formula for the 
Coulomb interaction energy, blue points include additional 
correction described in the text. Lines are least-squares fits 
with E(N) =a/N + E x . 
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Appendix A: Cohesive energy of MnO 

We have evaluated DMC energies of MnO crystal in 
three simulation cells containing 16, 24 and 32 atoms. 
The results obtained with the PBEOu, = o.3 orbitals are 
plotted in Fig. [5] Clearly, the averaging over 8 k points 
performed in all cases successfully suppresses differences 
in shapes of the individual simulation cells and the re- 
maining size dependence can be very well fitted with a 
function E(N) — a/N + E^. The estimate for E^ is 
-120.203 ±0.001 hartree/MnO. 

The a/N term can be substantially reduced by tech- 
niques described in Refs. [28| and HiJ One part of the 
proposed correction to the Ewald total energy reads 

^•■i/*T' (A1) 

The static structure factor S(k) entering the expression is 
defined as 5(k) = (^ol/SkP-kl^o) /N with standing for 
a Fourier component of the electron density. The integral 
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FIG. 6: (color online) Static structure factor S(k) in MnO 
(Bl AFM-II phase at V = 21.7 A 3 /MnO). Calculations were 
performed with the PBE0uj=o.3 orbitals. Data are combined 
from simulation cells containing 16, 24 and 32 atoms. Plot- 
ted are: VMC estimate (red squares), mixed DMC estimate 
(green empty circles) and extrapolated DMC estimate (blue 
filled circles). The lines are fits with S(k) = [l — exp(— ak 13 )] . 



in Eq. (|A1|) runs over a domain D centered around k = 
and having volume 87r 3 /f2, where f2 is volume of the sim- 
ulation cell. The structure factor S(k) is evaluated along 
the DMC simulation at a discrete set of points and then 
extrapolated towards k — 0. We show this extrapolation 
for the present case in Fig. [51 where we compare mixed 
and extrapolated DMC estimates, Eqs. ([2]) and ((3|). The 
small-fc behavior of S(k) is found to be ~ fc 17 for the 
mixed estimate and ~ fc 19 for the extrapolated one. This 
variance originates in a limited quality of our trial wave 
functions at large distances between electrons — the Jas- 
trow factor we use is restricted to zero for inter-electron 
separations larger than the Wigner-Seitz radius, and 
hence the correct long-range asymptotics cannot be fully 
captured. In the end, this deficiency is irrelevant, since 
the prediction of the extrapolated estimate is sufficiently 
close to the exact asymptotics ~ k 2 that translates to 
AEs(k) ~ l/jV . 28 i 29 The total energies corrected accord- 
ing to Eq. (jAip are shown in Fig. [5] together with the 
pure Ewald data. The correction reduces the finite-size 
errors by 75%. 

For calculations of individual atoms we utilize the same 
form of the trial wave function as we did in solids in order 
to stay within the same level of theoretical description. 



TABLE II: Fixed- node DMC total energy of isolated Mn atom 
calculated using several sets of single-particle orbitals. 

single-particle orbitals E (hartree) 

restricted open-shell HF -104.0133 ± 0.0006 
unrestricted HF -104.0185 ± 0.0005 

unrestricted PBE0„ =o .3 -104.0192 ± 0.0002 
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Single-particle orbitals for atomic calculations were ob- 
tained using Gamess code4^ Fixed-node DMC energies 
for the manganese atom corresponding to several choices 
of orbitals are listed in Tab. HU The difference between 
Hartree-Fock and PBEOm^o.3 is negligible for our pur- 
poses, since it is smaller than the statistical error bar 



achieved for the energy of the bulk crystal. 

Variation of single-particle orbitals in the trial wave 
function for the oxygen atom does not lead to any ap- 
preciable differences of DMC energies. To evaluate the 
cohesion of MnO, we use the result quoted in Ref. HH, 
E = -15.8421 ± 0.0002 hartree. 



* On leave from Institute of Physics, Academy of Sciences 
of the Czech Republic, Na Slovance 2, CZ-18221 Praha 8, 
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